library(psych)
library(tidyverse)
library(haven)
library(labelled)
library(car)
library(stargazer)
library(lfe)
library(survey)
library(readxl)
library(corrplot)
library(xtable)
library(glue)

setwd(paste0(dirname(rstudioapi::getActiveDocumentContext()$path)))

# functions needed for certain tests ####


difftest_felm <- function(x1, x2, model){
  diffest <- summary(model)$coef[x1,"Estimate"]-summary(model)$coef[x2,"Estimate"]
  vardiff <- (summary(model)$coef[x1,"Cluster s.e."]^2 + 
                summary(model)$coef[x2,"Cluster s.e."]^2) - (2*(vcov(model)[x1, x2])) 
  # variance of x1 + variance of x2 - 2*covariance of x1 and x2
  diffse <- sqrt(vardiff)
  tdiff <- (diffest)/(diffse)
  ptdiff <- 2*(1-pt(abs(tdiff), model$df, lower.tail=T))
  upr <- diffest + qt(.975, df = model$df)*diffse # will usually be very close to 1.96
  lwr <- diffest + qt(.025, df = model$df)*diffse
  df <- model$df
  return(c(est=round(diffest, digits =4), 
           diffse=round(diffse, digits = 4),
           t=round(tdiff, digits = 3), 
           p=round(ptdiff, digits = 4), 
           lwr=round(lwr, digits = 4), 
           upr=round(upr, digits = 4),
           df = df))
}

coef_star_test = function(mod,mod_std,var1,var2){
  dif = ifelse(coef(mod)[[var1]] >= coef(mod)[[var2]],"Yes","No")
  if(class(mod) == "felm"){
    test = difftest_felm(var1, var2, mod_std)
    star = ifelse(test[["p"]] <= 0.05, "*","")
  } else if(class(mod)=="lm"){
    test = linearHypothesis(mod_std,paste(var1,var2,sep="="))
    star = ifelse(test[2,ncol(test)]< 0.05, "*", "")
  } else{
    break
  }
  return(paste0(dif,star))
}

run_weighted_model <- function(formula, weights, data) {
  # first let's get the weights
  formula_vars <- all.vars(formula)
  data$complete_case <- as.integer(apply(data[formula_vars], 1, function(x) all(!is.na(x))))
  w_var = data[[weights]]
  full_weights <- 1 / table(w_var[data$complete_case==1])[as.character(w_var)]
  full_weights[is.na(full_weights)]=0
  # now let's run the regression
  mod = felm(formula, 
             weights = full_weights,
             data = data )
  return(mod)
}



# Load in original datasets ####

newspaper_dat = read_xlsx("welfare newspaper hits 90s v2.xlsx") # for figure 2
newspaper_dat_trans_raw = read_xlsx("trans newspaper hits 10s.xlsx") # for figure 3
trans_searches_raw = read.csv("transgender meaning google trends.csv", # for figure 4
                              col.names = c("transgender_meaning"))

# Read ANES datasets
cross_sectional_anes <- readRDS("cross_sectional_anes.rds") # for table 1
cross_sectional_anes_std <- readRDS("cross_sectional_anes_std.rds") # for table 1
two_year_panels_anes <- readRDS("two_year_panels_anes.rds") # for table 2
four_year_panels_anes <- readRDS("four_year_panels_anes.rds") # for table 2
two_year_panels_anes_std <- readRDS("two_year_panels_anes_std.rds") # for table 2
four_year_panels_anes_std <- readRDS("four_year_panels_anes_std.rds") # for table 2
anes_9297_sub <- readRDS("anes_9297_sub.rds") # for table 3
anes_9297_sub_std <- readRDS("anes_9297_sub_std.rds") # for table 3
anes_1620_sub <- readRDS("anes_1620_sub.rds") # for table 4
anes_1620_sub_std <- readRDS("anes_1620_sub_std.rds") # for table 4

# Read GSS datasets
cross_sectional_gss <- readRDS("cross_sectional_gss.rds") # for table 1
cross_sectional_gss_std <- readRDS("cross_sectional_gss_std.rds") # for table 1
two_year_panels_gss <- readRDS("two_year_panels_gss.rds") # for table 2
two_year_panels_gss_std <- readRDS("two_year_panels_gss_std.rds") # for table 2
four_year_panels_gss <- readRDS("four_year_panels_gss.rds") # for table 2
four_year_panels_gss_std <- readRDS("four_year_panels_gss_std.rds") # for table 2

# Results ####

setwd("../output")

## Table 1: Do values correlate with policy attitudes within-voter? ####

# ANES, trad 

summary(trad_crosssec_mod <- run_weighted_model(
  weights = "year",
  formula = trad_policies_2 ~ trad_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year  | 0 | VCASEID,
  
  data = cross_sectional_anes %>% 
    rename(trad_policies_2 = trad_policies,
           trad_values_1 = trad_values,
           trad_groups_1 = trad_groups,
           ideo_1 = ideo,
           pid_1 = pid,
           race_1 = race,
           educ_1 = educ,
           sex_1 = sex,
           inc_1 = inc)))


summary(trad_crosssec_mod_std <- run_weighted_model(
  weights = "year",
  formula = trad_policies_2 ~ trad_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year  | 0 | VCASEID,
  
  data = cross_sectional_anes_std %>% 
    rename(trad_policies_2 = trad_policies,
           trad_values_1 = trad_values,
           trad_groups_1 = trad_groups,
           ideo_1 = ideo,
           pid_1 = pid,
           race_1 = race,
           educ_1 = educ,
           sex_1 = sex,
           inc_1 = inc)))

summary(trad_2way_mod <- run_weighted_model(
  weights = "year",
  formula = trad_policies_2 ~ trad_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + VCASEID | 0 | VCASEID,
  
  data = cross_sectional_anes %>% 
    rename(trad_policies_2 = trad_policies,
           trad_values_1 = trad_values,
           trad_groups_1 = trad_groups,
           ideo_1 = ideo,
           pid_1 = pid,
           race_1 = race,
           educ_1 = educ,
           sex_1 = sex,
           inc_1 = inc)))


summary(trad_2way_mod_std <- run_weighted_model(
  weights = "year",
  formula = trad_policies_2 ~ trad_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + VCASEID | 0 | VCASEID,
  
  data = cross_sectional_anes_std %>% 
    rename(trad_policies_2 = trad_policies,
           trad_values_1 = trad_values,
           trad_groups_1 = trad_groups,
           ideo_1 = ideo,
           pid_1 = pid,
           race_1 = race,
           educ_1 = educ,
           sex_1 = sex,
           inc_1 = inc)))
linearHypothesis(trad_2way_mod_std,"trad_values_1=ideo_1")
linearHypothesis(trad_2way_mod_std,"trad_values_1=pid_1")

# ANES egal

summary(egal_crosssec_mod <- run_weighted_model(egal_policies_2 ~ egal_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year  | 0 | VCASEID, 
                                                weights="year",
                                                data = cross_sectional_anes %>% mutate(ideo = 1-ideo, pid = 1 - pid) %>% 
                                                  rename(egal_values_1 = egal_values,
                                                         egal_policies_2 = egal_policies,
                                                         egal_groups_1 = egal_groups,
                                                         ideo_1 = ideo,
                                                         pid_1 = pid,
                                                         race_1 = race,
                                                         educ_1 = educ,
                                                         sex_1 = sex,
                                                         inc_1 = inc)))

summary(egal_crosssec_mod_std <- run_weighted_model(egal_policies_2 ~ egal_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year | 0 | VCASEID, 
                                                    weights="year",
                                                    data = cross_sectional_anes_std %>% mutate(ideo = 1-ideo, pid = 1 - pid) %>% 
                                                      rename(egal_values_1 = egal_values,
                                                             egal_policies_2 = egal_policies,
                                                             egal_groups_1 = egal_groups,
                                                             ideo_1 = ideo,
                                                             pid_1 = pid,
                                                             race_1 = race,
                                                             educ_1 = educ,
                                                             sex_1 = sex,
                                                             inc_1 = inc)))
linearHypothesis(egal_crosssec_mod_std,"egal_values_1=ideo_1")
linearHypothesis(egal_crosssec_mod_std,"egal_values_1=pid_1")

summary(egal_2way_mod <- run_weighted_model(egal_policies_2 ~ egal_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + VCASEID | 0 | VCASEID, 
                                            weights="year",
                                            data = cross_sectional_anes %>% mutate(ideo = 1-ideo, pid = 1 - pid) %>% 
                                              rename(egal_values_1 = egal_values,
                                                     egal_policies_2 = egal_policies,
                                                     egal_groups_1 = egal_groups,
                                                     ideo_1 = ideo,
                                                     pid_1 = pid,
                                                     race_1 = race,
                                                     educ_1 = educ,
                                                     sex_1 = sex,
                                                     inc_1 = inc)))


summary(egal_2way_mod_std <- run_weighted_model(egal_policies_2 ~ egal_values_1 + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + VCASEID | 0 | VCASEID, 
                                                weights="year",
                                                data = cross_sectional_anes_std %>% mutate(ideo = 1-ideo, pid = 1 - pid) %>% 
                                                  rename(egal_values_1 = egal_values,
                                                         egal_policies_2 = egal_policies,
                                                         egal_groups_1 = egal_groups,
                                                         ideo_1 = ideo,
                                                         pid_1 = pid,
                                                         race_1 = race,
                                                         educ_1 = educ,
                                                         sex_1 = sex,
                                                         inc_1 = inc)))
linearHypothesis(egal_2way_mod_std,"egal_values_1=ideo_1")
linearHypothesis(egal_2way_mod_std,"egal_values_1=pid_1")

# GSS, trad

summary(trad_crosssec_mod_gss <- run_weighted_model(trad_policies2_gss ~ trad_values_1_gss + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year  | 0 | idnum, 
                                                    weights="year",data = cross_sectional_gss %>% 
                                                      rename(trad_policies2_gss = trad_policies,
                                                             trad_values_1_gss = trad_values,
                                                             ideo_1 = ideo,
                                                             pid_1 = pid,
                                                             race_1 = race,
                                                             educ_1 = educ,
                                                             sex_1 = sex,
                                                             inc_1 = inc)))



summary(trad_crosssec_mod_gss_std <- run_weighted_model(trad_policies2_gss ~ trad_values_1_gss + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year  | 0 | idnum, 
                                                        weights="year",data = cross_sectional_gss_std %>% 
                                                          rename(trad_policies2_gss = trad_policies,
                                                                 trad_values_1_gss = trad_values,
                                                                 ideo_1 = ideo,
                                                                 pid_1 = pid,
                                                                 race_1 = race,
                                                                 educ_1 = educ,
                                                                 sex_1 = sex,
                                                                 inc_1 = inc)))


summary(trad_2way_mod_gss <- run_weighted_model(trad_policies2_gss ~ trad_values_1_gss + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + idnum | 0 | idnum, 
                                                weights="year",data = cross_sectional_gss %>% 
                                                  rename(trad_policies2_gss = trad_policies,
                                                         trad_values_1_gss = trad_values,
                                                         ideo_1 = ideo,
                                                         pid_1 = pid,
                                                         race_1 = race,
                                                         educ_1 = educ,
                                                         sex_1 = sex,
                                                         inc_1 = inc)))

summary(trad_2way_mod_gss_std <- run_weighted_model(trad_policies2_gss ~ trad_values_1_gss + ideo_1 + pid_1 + as.factor(race_1) + as.factor(educ_1) + sex_1 + inc_1 | year + idnum | 0 | idnum, 
                                                    weights="year",data = cross_sectional_gss_std %>% 
                                                      rename(trad_policies2_gss = trad_policies,
                                                             trad_values_1_gss = trad_values,
                                                             ideo_1 = ideo,
                                                             pid_1 = pid,
                                                             race_1 = race,
                                                             educ_1 = educ,
                                                             sex_1 = sex,
                                                             inc_1 = inc)))

linearHypothesis(trad_2way_mod_gss_std,"trad_values_1_gss=ideo_1")
linearHypothesis(trad_2way_mod_gss_std,"trad_values_1_gss=pid_1")

within_mod_sg = stargazer(trad_crosssec_mod,trad_2way_mod,
                          egal_crosssec_mod,egal_2way_mod,
                          trad_crosssec_mod_gss,trad_2way_mod_gss,
                          title = "Cross-sectional relationships between values and policy views",
                          label="within_mod_sg",
                          single.row = F,
                          omit = c("prediction","race","sex","inc","educ") ,
                          
                          covariate.labels = c("Trad. Values (ANES)",
                                               "Egal. Values (ANES)",
                                               "Trad. Values (GSS)",
                                               "Ideology",
                                               "Partisanship"),
                          
                          
                          column.sep.width = "-8pt",
                          no.space=T,
                          table.placement = "!hbtp",
                          dep.var.labels = c("\\multirow{2}{4cm}{\\centering Cultural Orthodoxy Policy Index (ANES)}",
                                             "\\multirow{2}{4cm}{\\centering Opportunity Policy Index (ANES)}",
                                             "\\multirow{2}{4cm}{\\centering Cultural Orthodoxy Policy Index (GSS)}"),
                          add.lines = list(c("Respondent Fixed Effects",rep(c("\\xmark","\\checkmark"),3) ),
                                           c("Values > Ideology",
                                             coef_star_test(trad_crosssec_mod,trad_crosssec_mod_std,"trad_values_1","ideo_1"),
                                             coef_star_test(trad_2way_mod,trad_2way_mod_std,"trad_values_1","ideo_1"),
                                             coef_star_test(egal_crosssec_mod,egal_crosssec_mod_std,"egal_values_1","ideo_1"),
                                             coef_star_test(egal_2way_mod,egal_2way_mod_std,"egal_values_1","ideo_1"),
                                             coef_star_test(trad_crosssec_mod_gss,trad_crosssec_mod_gss_std,"trad_values_1_gss","ideo_1"),
                                             coef_star_test(trad_2way_mod_gss,trad_2way_mod_gss_std,"trad_values_1_gss","ideo_1")),
                                           c("Values > Partisanship",
                                             coef_star_test(trad_crosssec_mod,trad_crosssec_mod_std,"trad_values_1","pid_1"),
                                             coef_star_test(trad_2way_mod,trad_2way_mod_std,"trad_values_1","pid_1"),
                                             coef_star_test(egal_crosssec_mod,egal_crosssec_mod_std,"egal_values_1","pid_1"),
                                             coef_star_test(egal_2way_mod,egal_2way_mod_std,"egal_values_1","pid_1"),
                                             coef_star_test(trad_crosssec_mod_gss,trad_crosssec_mod_gss_std,"trad_values_1_gss","pid_1"),
                                             coef_star_test(trad_2way_mod_gss,trad_2way_mod_gss_std,"trad_values_1_gss","pid_1"))),
                                                table.layout ="-ldc-#-t-a-s-n",
                          keep.stat = c("n","adj.rsq"),
                          notes.align = "l",
                          star.cutoffs = 0.05,star.char = "*",
                          notes.append = FALSE,
                          notes.label = "",
                          notes = "\\parbox[t]{6.5in}{\\footnotesize \\textit{Note}: Data weighted to give each survey wave equal weight. All models include survey wave fixed effects and demographic controls. Standard errors clustered at the respondent level. 
                          Ideology and partisanship variables are coded such that higher values are more conservative/Republican for cultural orthodoxy policy models and more liberal/Democratic for opportunity policy models. $^{*}$ indicates \\emph{p} $<$ 0.05 (two-tailed tests). }")
within_mod_sg[12] = paste0("\\\\",within_mod_sg[12])

cat(within_mod_sg, sep = '\n', file = "within_mod_sg.tex")

## Table 2: Do values drive changes in policy views? ####


summary(trad_2year_mod <- run_weighted_model(trad_policies_2 ~ trad_values_1 + trad_policies_1 + ideo_1 + pid_1  | year1 | 0 | VCASEID, 
                                             weights = "year1", data = two_year_panels_anes ))
summary(trad_2year_mod_std <- run_weighted_model(weights = "year1", formula = trad_policies_2 ~ trad_values_1 + trad_policies_1 + ideo_1 + pid_1  | year1 | 0 | VCASEID, 
                                                 data = two_year_panels_anes_std))

difftest_felm("trad_values_1", "ideo_1", trad_2year_mod_std)
difftest_felm("trad_values_1", "pid_1", trad_2year_mod_std)

summary(trad_4year_mod <- run_weighted_model(weights = "year1", formula = trad_policies_2 ~ trad_values_1 + trad_policies_1 + ideo_1 + pid_1  | year1 | 0 | VCASEID, data = four_year_panels_anes))
summary(trad_4year_mod_std <- run_weighted_model(weights = "year1", formula = trad_policies_2 ~ trad_values_1 + trad_policies_1 + ideo_1 + pid_1  | year1 | 0 | VCASEID, data = four_year_panels_anes_std))

difftest_felm("trad_values_1", "ideo_1", trad_4year_mod_std)
difftest_felm("trad_values_1", "pid_1", trad_4year_mod_std)



summary(egal_2year_mod <- run_weighted_model(weights = "year1", formula = egal_policies_2 ~ egal_values_1 +egal_policies_1 +  ideo_1 + pid_1  | year1 | 0 | VCASEID, 
                                             data = two_year_panels_anes %>% mutate(ideo_1 = 1-ideo_1, pid_1 = 1 - pid_1)))
summary(egal_2year_mod_std <- run_weighted_model(weights = "year1", formula = egal_policies_2 ~ egal_values_1 +egal_policies_1 +  ideo_1 + pid_1  | year1 | 0 | VCASEID, 
                                                 data = two_year_panels_anes_std %>% mutate(ideo_1 = 1-ideo_1, pid_1 = 1 - pid_1)))

difftest_felm("egal_values_1", "ideo_1", egal_2year_mod_std)
difftest_felm("egal_values_1", "pid_1", egal_2year_mod_std)

summary(egal_4year_mod <- run_weighted_model(weights = "year1", formula = egal_policies_2 ~ egal_values_1 +egal_policies_1 +  ideo_1 + pid_1 | year1 | 0 | VCASEID, 
                                             data = four_year_panels_anes %>% mutate(ideo_1 = 1-ideo_1, pid_1 = 1 - pid_1)))
summary(egal_4year_mod_std <- run_weighted_model(weights = "year1", formula = egal_policies_2 ~ egal_values_1 +egal_policies_1 +  ideo_1 + pid_1 | year1 | 0 | VCASEID,
                                                 data = four_year_panels_anes_std %>% mutate(ideo_1 = 1-ideo_1, pid_1 = 1 - pid_1)))

difftest_felm("egal_values_1", "ideo_1", egal_4year_mod_std)
difftest_felm("egal_values_1", "pid_1", egal_4year_mod_std)


summary(trad_2year_mod_gss <- run_weighted_model(weights = "year1", formula = trad_policies_2_gss ~ trad_values_1_gss +trad_policies_1_gss +  ideo_1 + pid_1 | year1 | 0 | idnum, 
                                                 data = two_year_panels_gss %>% 
                                                   rename(trad_policies_2_gss = trad_policies_2,
                                                          trad_policies_1_gss = trad_policies_1,
                                                          trad_values_1_gss = trad_values_1
                                                   )))
summary(trad_2year_mod_std_gss <- run_weighted_model(weights = "year1", formula = trad_policies_2 ~ trad_policies_1 + trad_values_1_gss + ideo_1 + pid_1 | year1 | 0 | idnum, 
                                                     data = two_year_panels_gss_std %>% mutate(trad_values_1_gss = trad_values_1)))

difftest_felm("trad_values_1_gss", "ideo_1", trad_2year_mod_std_gss)
difftest_felm("trad_values_1_gss", "pid_1", trad_2year_mod_std_gss)

summary(trad_4year_mod_gss <-run_weighted_model(weights = "year1", formula = trad_policies_2_gss ~ trad_values_1_gss + trad_policies_1_gss +  ideo_1 + pid_1 | year1 | 0 | idnum, 
                                                data = four_year_panels_gss %>% 
                                                  rename(trad_policies_2_gss = trad_policies_2,
                                                         trad_policies_1_gss = trad_policies_1,
                                                         trad_values_1_gss = trad_values_1
                                                  )))
summary(trad_4year_mod_std_gss <- run_weighted_model(weights = "year1", formula = trad_policies_2 ~ trad_policies_1 + trad_values_1_gss + ideo_1 + pid_1 | year1 | 0 | idnum, data = four_year_panels_gss_std %>% rename(trad_values_1_gss = trad_values_1)))

difftest_felm("trad_values_1_gss", "ideo_1", trad_4year_mod_std_gss)
difftest_felm("trad_values_1_gss", "pid_1", trad_4year_mod_std_gss)


shift_mod_sg = stargazer(trad_2year_mod,trad_4year_mod,egal_2year_mod,egal_4year_mod,
                         trad_2year_mod_gss,trad_4year_mod_gss,
                         title = "Relationship between values and changes in policy views",
                         label="shift_mod_sg",
                         single.row = F,
                         order = c(1,3,5,7,8,2,4,6),
                         covariate.labels = c("Traditionalist values (ANES)",
                                              "Egalitarian values (ANES)",
                                              "Traditionalist values (GSS)",
                                              "Ideology",
                                              "Partisanship",
                                              "Cultural Orth. Policy Index (ANES)",
                                              "Opportunity Policy Index (ANES)",
                                              "Cultural Orth. Policy Index (GSS)"),
                         
                         column.sep.width = "-15pt",
                         no.space=T,
                         table.placement = "!hbtp",
                         omit = "prediction",
                         add.lines = list(c("Values > Ideology",
                                            coef_star_test(trad_2year_mod,trad_2year_mod_std,"trad_values_1","ideo_1"),
                                            coef_star_test(trad_4year_mod,trad_4year_mod_std,"trad_values_1","ideo_1"),
                                            coef_star_test(egal_2year_mod,egal_2year_mod_std,"egal_values_1","ideo_1"),
                                            coef_star_test(egal_4year_mod,egal_4year_mod_std,"egal_values_1","ideo_1"),
                                            coef_star_test(trad_2year_mod_gss,trad_2year_mod_std_gss,"trad_values_1_gss","ideo_1"),
                                            coef_star_test(trad_2year_mod_gss,trad_4year_mod_std_gss,"trad_values_1_gss","ideo_1")),
                                          c("Values > Partisanship",
                                            coef_star_test(trad_2year_mod,trad_2year_mod_std,"trad_values_1","pid_1"),
                                            coef_star_test(trad_4year_mod,trad_4year_mod_std,"trad_values_1","pid_1"),
                                            coef_star_test(egal_2year_mod,egal_2year_mod_std,"egal_values_1","pid_1"),
                                            coef_star_test(egal_4year_mod,egal_4year_mod_std,"egal_values_1","pid_1"),
                                            coef_star_test(trad_2year_mod_gss,trad_2year_mod_std_gss,"trad_values_1_gss","pid_1"),
                                            coef_star_test(trad_2year_mod_gss,trad_4year_mod_std_gss,"trad_values_1_gss","pid_1"))),
                         
                         
                         dep.var.labels = c("\\multirow{2}{4cm}{\\centering Cultural Orthodoxy Policy Index (ANES)}",
                                            "\\multirow{2}{4cm}{\\centering Opportunity Policy Index (ANES)}", 
                                            "\\multirow{2}{4cm}{\\centering Cultural Orthodoxy Policy Index (GSS)}"),
                         column.labels = rep(c("+2 yrs",
                                               "+4 yrs"),4),
                         table.layout ="-ldc-#-t-a-s-n",
                         keep.stat = c("n","adj.rsq"),
                         notes.align = "l",
                         star.cutoffs = 0.05,star.char = "*",
                         notes.append = FALSE,
                         notes.label = "",
                         notes = "\\parbox[t]{6.7in}{\\footnotesize \\textit{Note}: ANES data pooled from three two-year panels (92-94, 94-96, 00-02) and three four-year panels (92-96, 00-04, 16-20). GSS data pooled from five two-year panels (06-08, 08-10, 10-12, 12-14, 18-20) and three four-year panels (06-10, 10-14, 16-20). 
                         Data weighted to give each survey wave equal weight. All models include survey wave fixed effects. Standard errors clustered at the respondent level. 
                          Ideology and partisanship variables are coded such that higher values are more conservative/Republican for cultural orthodoxy policy models and more liberal/Democratic for opportunity policy models. $^{*}$ indicates \\emph{p} $<$ 0.05 (two-tailed tests).}")
shift_mod_sg[12] = paste0("\\\\",shift_mod_sg[12])

cat(shift_mod_sg, sep = '\n', file = "shift_mod_sg.tex")

## Figure 2: Newspaper plot, welfare reform ####

newspaper_dat_adj = newspaper_dat %>%
  mutate_at(vars(congress:welfare_immigrant_eligib), function(x) x / newspaper_dat$monday) %>%
  arrange(year) 

newspaper_dat_adj2 = newspaper_dat_adj %>% 
  sweep(2, unlist(newspaper_dat_adj[1, ]), `/`) %>%
  mutate(year = year * 1990)

newspaper_dat_adj2 %>% 
  select(-monday) %>%
  pivot_longer(-year,values_to="val",names_to="name") %>%
  mutate(name = case_when(
    name == "family_cap" ~ '"family cap"',
    name == "welfare_year_limit" ~ 'welfare + "year limit"',
    name == "welfare_work_requirement" ~ 'welfare + "work requirement"',
    name == "welfare_immigrant_eligib" ~ 'welfare + immigrant + eligib*',
    T ~ name)) %>% 
  mutate(Term = fct_relevel(as.factor(name),"congress","welfare",'welfare + "year limit"',
                            'welfare + "work requirement"','welfare + immigrant + eligib*',
                            '"family cap"')) %>% 
  ggplot(aes(x=year,y=val,color=Term))+
  geom_line(size=1)+
  theme_minimal()+
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA))+
  scale_x_continuous(breaks = 1990:1999)+
  scale_y_continuous(trans='log2',breaks=c(0.5,1,2,4,8,16,32),
                     labels = function(x) str_c(x, 'x'))+
  labs(x="Year",y=str_wrap("Newspaper Coverage Relative to 1990 (Log Scale)",35))

ggsave("welfare_term_usage.png",width=7,height=5,dpi=600)

## Table 3: Welfare policy attitudes ####

summary(egal_92_mod_welfare <- lm(welfare_policies_96 ~ egal_values_92 + ideo__92 +  pid__92 + welfare_policies_92 , 
                                  data = anes_9297_sub %>% 
                                    mutate_at(vars(starts_with("ideo")), function(x) 1 - x ) %>% 
                                    mutate_at(vars(starts_with("pid")), function(x) 1 - x )))

summary(egal_92_mod_welfare_std <- lm(welfare_policies_96 ~ egal_values_92 + ideo__92 +  pid__92 + welfare_policies_92 , 
                                      data = anes_9297_sub_std %>% 
                                        mutate_at(vars(starts_with("ideo")), function(x) 1 - x ) %>% 
                                        mutate_at(vars(starts_with("pid")), function(x) 1 - x )))
linearHypothesis(egal_92_mod_welfare_std, "egal_values_92 = pid__92")

summary(egal_94_mod_welfare <- lm(welfare_policies_96 ~ egal_values_94 + ideo__94 + pid__94 +  welfare_policies_94 ,
                                  data = anes_9297_sub %>% 
                                    mutate_at(vars(starts_with("ideo")), function(x) 1 - x ) %>% 
                                    mutate_at(vars(starts_with("pid")), function(x) 1 - x )))
summary(egal_94_mod_welfare_std <- lm(welfare_policies_96 ~ egal_values_94 + ideo__94 + pid__94 +  welfare_policies_94 ,
                                      data = anes_9297_sub_std %>% 
                                        mutate_at(vars(starts_with("ideo")), function(x) 1 - x ) %>% 
                                        mutate_at(vars(starts_with("pid")), function(x) 1 - x )))
linearHypothesis(egal_94_mod_welfare_std, "egal_values_94 = ideo__94")
linearHypothesis(egal_94_mod_welfare_std, "egal_values_94 = pid__94")


welfare_sg = stargazer(egal_92_mod_welfare,egal_94_mod_welfare,
                       title = "Relationship between values and future views on welfare reform",
                       label="welfare_sg",
                       single.row = T,
                       covariate.labels = c("Egalitarian Values (1992)",
                                            "Ideology (1992)",
                                            "Partisanship (1992)",
                                            "Welfare Policy Index (1992)",
                                            "Egalitarian Values (1994)",
                                            "Ideology (1994)",
                                            "Partisanship (1994)",
                                            "Welfare Policy Index (1994)"),
                       column.sep.width = "-5pt",
                       no.space=T,
                       table.placement = "!hbtp",
                       omit = "Constant",
                       add.lines = list(c("Egal. Values > Ideology",
                                          coef_star_test(egal_92_mod_welfare,egal_92_mod_welfare_std,"egal_values_92","ideo__92"),
                                          coef_star_test(egal_94_mod_welfare,egal_94_mod_welfare_std,"egal_values_94","ideo__94")),
                                        c("Egal. Values > Partisanship",
                                          coef_star_test(egal_92_mod_welfare,egal_92_mod_welfare_std,"egal_values_92","pid__92"),
                                          coef_star_test(egal_94_mod_welfare,egal_94_mod_welfare_std,"egal_values_94","pid__94"))),
                       dep.var.labels = c("Welfare Policy Index (1996)"),
                       table.layout ="-ld-#-t-as-n",
                       keep.stat = c("n","adj.rsq"),
                       notes.align = "l",
                       star.cutoffs = 0.05,star.char = "*",
                       notes.append = FALSE,
                       notes.label = "",
                       notes = "\\parbox[t]{5in}{\\footnotesize \\textit{Note}: $^{*}$ indicates \\emph{p} $<$ 0.05 (two-tailed tests). Ideology and partisanship variables are coded such that higher values are more liberal/Democratic. }")

cat(welfare_sg, sep = '\n', file = "welfare_sg.tex")

## Figure 3: Newspaper frequency, transgender-related terms ####

newspaper_dat_trans <- newspaper_dat_trans_raw %>%
  mutate(across(congress:transgender_military, ~ . / cur_data()[["monday"]])) %>%
  arrange(year)

newspaper_dat_trans2 = newspaper_dat_trans %>% 
  sweep(2, unlist(newspaper_dat_trans[1, ]), `/`) %>%
  mutate(year = year * 2010)

newspaper_dat_trans2 %>% 
  select(-monday) %>%
  pivot_longer(-year,values_to="val",names_to="name") %>%
  drop_na() %>% 
  mutate(name = case_when(
    name == "gender_identity" ~ '"gender identity"',
    name == "transgender_bathroom" ~ 'transgender + bathroom',
    name == "transgender_military" ~ 'transgender + military',
    name == "transgender_rights" ~ 'transgender + rights',
    T ~ name)) %>%
  mutate(Term = fct_relevel(as.factor(name),"congress","transgender",'transgender + rights',
                            'transgender + bathroom','transgender + military',
                            '"gender identity"')) %>%
  ggplot(aes(x=year,y=val,color=Term))+
  geom_line(size=1)+
  theme_minimal()+
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA))+
  scale_x_continuous(breaks = seq(2010,2022,2))+
  scale_y_continuous(trans='log2',breaks=c(0.5,1,2,4,8,16,32,64),
                     labels = function(x) str_c(x, 'x'))+
  labs(x="Year",y=str_wrap("Newspaper Coverage Relative to 2010 (Log Scale)",30),
       color = "")

ggsave("transgender_term_usage.png",width=7,height=5, dpi = 600)

## Figure 4: Google trends, transgender meaning ####

trans_searches = trans_searches_raw %>% 
  rownames_to_column("date") %>% 
  slice(-1) %>%
  mutate(date = as.Date(paste(date, "-01", sep=""), format="%Y-%m-%d"),
         transgender_meaning= as.integer(transgender_meaning))  %>%
  filter(date >= as.Date("2010-01-01") & date <= as.Date("2022-12-01"))


ggplot(trans_searches,
       aes(x=date,y=transgender_meaning)) + 
  geom_line(size=1)+
  theme_minimal()+
  theme(legend.position = "bottom",
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA))+ 
  scale_x_date(date_breaks = "2 years", date_labels = "%Y",
               expand = c(0, 0))+
  
  labs(x="Year",y=str_wrap('Google Trends Searches for "transgender" + "meaning" (indexed to max=100)',30))

ggsave("transgender_term_usage2.png",width=7,height=5,dpi=600)


## Table 4: Transgender policies ####

summary(trans_1620_mod2 <- lm(trad_transpolicies_20 ~ trad_values_16 + egal_values_16 + ideo__16 + pid__16 + trad_transpolicies_16, 
                              data = anes_1620_sub %>% mutate(egal_values_16 = 1 - egal_values_16 )))
summary(trans_1620_mod2_std <- lm(trad_transpolicies_20 ~ trad_values_16 +  egal_values_16+ ideo__16 + pid__16 + trad_transpolicies_16, 
                                  data = anes_1620_sub_std %>% mutate(egal_values_16 = 1 - egal_values_16 )))
linearHypothesis(trans_1620_mod2_std,"trad_values_16 = ideo__16")
linearHypothesis(trans_1620_mod2_std,"trad_values_16 = pid__16")
linearHypothesis(trans_1620_mod2_std,"egal_values_16 = ideo__16")
linearHypothesis(trans_1620_mod2_std,"egal_values_16 = pid__16")


trans_sg = stargazer(trans_1620_mod2,
                     title = "Relationship between values and future views on transgender policy",
                     label="trans_sg",
                     single.row = T,
                     covariate.labels = c("Traditional Values (2016)",
                                          "Egalitarian Values (2016, Reversed)",
                                          "Ideology (2016)",
                                          "Partisanship (2016)",
                                          "Trans. Bathroom (2016)"),
                     column.sep.width = "0pt",
                     no.space=T,
                     table.placement = "!hbtp",
                     omit = "Constant",
                     dep.var.labels = c("Trans. Policy Index (2020)"),
                     add.lines = list(c("Trad. Values > Ideology",
                                        coef_star_test(trans_1620_mod2,trans_1620_mod2_std,"trad_values_16","ideo__16")),
                                      c("Trad. Values > Partisanship",
                                        coef_star_test(trans_1620_mod2,trans_1620_mod2_std,"trad_values_16","pid__16")),
                                      c("Egal. Values > Ideology",
                                        coef_star_test(trans_1620_mod2,trans_1620_mod2_std,"egal_values_16","ideo__16")),
                                      c("Egal. Values > Partisanship",
                                        coef_star_test(trans_1620_mod2,trans_1620_mod2_std,"egal_values_16","pid__16"))),
                     table.layout ="-ld-#-t-as-n",
                     keep.stat = c("n","adj.rsq"),
                     notes.align = "l",
                     star.cutoffs = 0.05,star.char = "*",
                     notes.append = FALSE,
                     notes.label = "",
                     notes = "\\parbox[t]{5in}{\\footnotesize \\textit{Note}: Ideology and partisanship variables are coded such that higher values are more conservative/Republican. Egalitarianism variable is reverse coded such that higher values indicate lower degrees of egalitarianism. $^{*}$ indicates \\emph{p} $<$ 0.05 (two-tailed tests). }")


cat(trans_sg, sep = '\n', file = "trans_sg.tex")
